#include<iostream>
#include<cmath>
#include<vector>
using namespace std;
// Function to calculate the coefficients of the best fitting polynomial
// using the normal equations
vector<double> discreteLeastSquare(vector<double> x, vector<double> y) {
int n = x.size();
double sum_x = 0, sum_x2 = 0, sum_x3 = 0, sum_x4 = 0;
double sum_y = 0, sum_xy = 0, sum_x2y = 0;

// Calculate the sums
for (int i = 0; i < n; i++) {
    sum_x += x[i];
    sum_x2 += pow(x[i], 2);
    sum_x3 += pow(x[i], 3);
    sum_x4 += pow(x[i], 4);
    sum_y += y[i];
    sum_xy += x[i] * y[i];
    sum_x2y += pow(x[i], 2) * y[i];
}

// Create the G
double G[3][3] = {
{n, sum_x, sum_x2},
{sum_x, sum_x2, sum_x3},
{sum_x2, sum_x3, sum_x4}
};


// Create the vector
double C[3] = {sum_y, sum_xy, sum_x2y};

// Perform Gaussian elimination
for (int i = 0; i < 3; i++) {
double pivot = G[i][i];
for (int j = 0; j < 3; j++) {
G[i][j] /= pivot;
}
C[i] /= pivot;
for (int k = 0; k < 3; k++) {
if (k != i) {
double factor = G[k][i];
for (int j = 0; j < 3; j++) {
G[k][j] -= factor * G[i][j];
}
C[k] -= factor * C[i];
}
}
}

// Return the coefficients
vector<double> coefficients; 
for (int i = 0; i < 3; i++) {
coefficients.push_back(C[i]);
}
return coefficients;
}
int main() {
// Input arrays
vector<double> x = {0.0 ,0.5, 1.0, 1.5, 2.0, 2.5, 3.0,3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5,7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0};
vector<double> y = {2.9 ,2.7 ,4.8 ,5.3 ,7.1 ,7.6, 7.7,7.6, 9.4, 9.0, 9.6, 10.0, 10.2, 9.7,8.3, 8.4, 9.0, 8.3, 6.6, 6.7, 4.1};

// Calculate the coefficients
vector<double> coefficients = discreteLeastSquare(x, y);

// Print the coefficients
cout << "a0 = " << coefficients[0] << endl;
cout << "a1 = " << coefficients[1] << endl;
cout << "a2 = " << coefficients[2] << endl;

return 0;
}